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^1 It is well-known that the synchronization of diffusively-coupled systems on networks strongly 

■^ ■ depends on the network topology. In particular, the so-called algebraic connectivity fiN-i, or the 

(^ 

smallest non-zero eigenvalue of the discrete Laplacian operator plays a crucial role on synchroniza- 
tion, graph partitioning, and network robustness. In our study, synchronization is placed in the 
C/^ ' general context of networks-of-networks, where single network models are replaced by a more real 

Q ■ istic hierarchy of interdependent networks. The present work shows, analytically and numerically, 

how the algebraic connectivity experiences sharp transitions after the addition of sufficient links 
among interdependent networks. 
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I. INTRODUCTION 

In the past decades, there has been a significant advance in understanding the structure 
and function of complex networks [l|, |2] • Mathematical models of networks are now widely 
used to describe a broad range of complex systems, from networks of human contacts to in- 
teractions amongst proteins. In particular, synchronization as an emerging phenomenon of 
a population of dynamically interacting units has always fascinated humans. The scientific 
interest in synchronization of coupled oscillators can be traced back to the work by Chris- 
tiaan Huygens on "An odd kind sympathy", between coupled pendulum clocks |3|, where he 
noticed that two pendulum clocks mounted on the same frame will synchronize after some 
time. Synchronization phenomena and processes are ubiquitous in nature and play a vital 



role within various contexts in biology, chemistry, ecology, sociology, and technology 



Ml. 



a. 



To date, the problem 



extending topics as epidemic spread [5| and coupled oscillators 

of how the structural properties of a network influence the performance and stability of the 

fully synchronized states of the network have been extensively investigated and discussed. 



both numerically and theoretically |9l-ll2| 



There exist many different definitions for synchronization (most of which are related). 
However, in the present paper, we will employ a definition based on the following equations: 

jeNi 
where Si represents the (relative) deviation of the i—th component state from its equilibrium, 
Ni its neighbors, and Q the Laplacian matrix, as will be defined in section [TTl In words, we 
model each component by a differential equation, such that the equations of the whole 
system become coupled due to the linking of the components in the network. Assuming 
such a perspective, the synchronization of a network maps into the dynamics of ([1]). The 
paper has been mostly written with this type of application in mind. Nevertheless, results 
extend to all phenomena dominated by the Laplacian, such as diffusion delivery of any 
commodity on a network. 

It is well-known that the synchronization of diffusively-coupled systems on networks is 



crucially affected by the network topology |l3l-|l6|. However, current research methods 



focus almost exclusively on individual networks treated as isolated systems. In reality, an 
individual network is often a combined system of multiple networks with distinct topologies 



and functions. This motivates us to study the effect of interdependent topologies on the 
mutual synchronization of networks. Recently, effort has been directed to these complex 
systems composed of many interdependent networks, which seem to model complex systems 



better than single networks 17|,|l8|. For instance, a pathogen spreads on a network of human 



contacts supported by global and regional transportation networks; or in a power grid and 



a communication network, that are coupled together [15|, a power station depends on a 
communication node for control, while a communication node depends on a power station 
for electricity. Cascading failures on interdependent networks, where the failure of a node 
at one end of an interdependent link implies the failure of the node at the other end of the 



link, have been widely studied [18|, |19|. The latter studies show that results obtained in the 
context of a single isolated network can change dramatically once interactions with other 
networks are incorporated. 

In particular, we will focus on the the so-called algebraic connectivity of interdependent 
networks, which is defined as the second smallest eigenvalue fiN-i of the discrete Lapla- 
cian matrix. This eigenvalue plays an important role on, among others, synchronization 
dynamics, network robustness, consensus problems, flocking and swarming, belief propa- 
gation, synchronization of coupled oscillators, graph partitioning, distributed filtering in 



sensor networks 



13-16 



20I-I24 



-24j. In the present work, we interpret the algebraic connectivity 
as the inverse of a "proper time", since the deviations from equilibrium in ([1]) decay ex- 
ponentially with such scale. Larger values of fiN-i enable synchronization in both discrete 



and continuous-time systems, even in the presence of transmission delays 20|, |25[. From 
a graph theoretic perspective, we will show that the algebraic connectivity experiences a 
phase transition upon the addition of a sufficient number of links among two interdependent 
networks. In other words, system synchronizability does not experience any depletion when 
the operability of the control channel is softly reduced. 

This paper is structured as follows. Section [TTl introduces the necessary notation and 
exposes both the Laplacian matrix and the graph spectra. Sections IIIII and IIVI provide a 
mean-field approach for the algebraic connectivity, and exploit the perturbation theory of 
interdependent networks, respectively. Finally, numerical results are presented in Section IVl 
The latter will expose properties of regular, random, small-world, and scale-free networks. 
Conclusions are drawn in Section |Vll 



II. DEFINITIONS 



A. Graph Theory Basics 



A graph G is composed by a set of nodes interconnected by a set of links G(Af,C). 
Suppose one has two networks Gi = (A/i,£i) and G2 = (A/2, £2), each with a set of nodes 
(A/i, A/2) and a set of hnks (£1, £2) respectively. For simplicity, in this paper we only study 
the case where Gi and G2 are identical, i.e. Gi = G2, meaning that the i-th node of Gi 
is topologically equivalent to the i-th node in G2. In the following, we will suppose any 
dependence relation to be symmetric, i.e. all networks are undirected. 

The global system resulting from the connection of the two networks is a network G 
with A/i U A/2 nodes and Ci U £2 "intralinks" plus a number of "interlinks" C12 joining 
the two networks; that is A/" = A/i U J\f2 and £ = £1 U £2 U £12, thus (A/", C) = G = 
(AriUAr2,£iU£2U£i2). 

Let us denote Ni as the number of nodes in lA/^l, and Lj as the number of links in as 
|£j|, also A^ = A''i + N2, and L = Li + L2; let Ai and A2 be the adjacency matrices of the 
two networks Gi and G2, and A that of the whole system G, whose entries or elements are 
ttij = 1 if node i is connected to node j, otherwise aij = 0. When the two networks are 
disconnected (£12 = 0), the matrix A is defined as the N x N matrix: 



A 



Ai 
A2 

When an interaction is introduced (£12 7^ 0), the adjacency matrix acquires non-trivial 
off-block terms denoted by -B/j, defined as the A/j x Nj interconnection matrix representing 
the interlinks between Gi and 6*2- The interdependency matrix B is then 



B 



B 



12 



BI2 



When the two networks Gi and G2 are equal, the adjacency matrix of the total system 
can be written as: 



A + aB 



Ai aB 



12 



aB{^ A2 



(2) 



where a represents coupling strength of the interaction. If the type of relation inside the 
A and B networks is the same (i.e. if they represent infrastructures of the same type, such 
as for instance electric systems), one may study the properties of non- weighted adjacency 
matrices. 

Similarly to the adjacency matrix, one may introduce the Laplacian matrix Q = D — A\ 
where D is the diagonal matrix of the degrees, where the degree of the i-th node is 
di = ^ • aij. In the same vein, one may define the diagonal matrices: 



(D 



def 



1 a 



EAB 



12jiJ5 



def 



{D2)ii "'= Yjj{B2l)ij = Yjji.Bi2)ii'^ 

and the Laplacian Q of the total system G reads: 



Q = Qa + oiQb 



(3) 



Qi + aDi -aBi2 
-aBf^ Q2 + aD2 

where Qi = Q2 is the Laplacian matrix of Ai = A2, and Qb is the Laplacian only repre- 
senting the interlinks: 



Q 



B 



D-B 



-BI2 



-B12 
D2 



(4) 



B. Fiedler Partitioning 

Since Q is a real symmetric matrix, it has A^ real eigenvalues |26|, which we order non- 
decreasingly = hm < jJ'N-i < ■ ■ ■ < /Ui- The eigenvector xm~i corresponding to the first 
non-zero eigenvalue ^n-i provides a graph partition named after Fiedler, who derived the 
majority of its properties J23|, |27|. The A^ — 1 largest Laplacian eigenvectors and eigenvalues 
satisfy the following equations: 



Qx 



T 
X X 



T 
X U 



fix, 
1, 

0. 



(5) 



where u is the all ones vector, which is the Laplacian eigenvector belonging to fi^ = 0. The 
algebraic connectivity f^N-i is the smallest of the N — 1 eigenvalues satisfying the equations 



in ([5]). Equivalently, xn-i and jj,n-i optimize the quadratic form x Qx subject to two 
constraints: 



Hn-i = min x Qx. (6) 

x'^=l,x'^u=0 

Since we will only deal with the Fiedler eigenvector, we will simplify the notation of the 
eigenpair (/i7v_i,X7v_i) by simply writing (/i,x). 



III. EXACT RESULTS FOR MEAN-FIELD THEORY 



A. Diagonal interlinking 



Let us start with the case of two exactly identical networks connected by C12 correspond- 
ing interlinks. The mean-field approach to such a system consists in studying a graph of two 
identical networks interacting via A'^i weighted connections among all corresponding nodes. 
The weight of each link, represented by a = ■^, equals the fraction of nodes linked to their 
corresponding in the exact network. In other words B12 = I, such that the synchronization 
interdependence is modulated by the parameter a: 



Q 



B 



I -I 
-I I 



(7) 



and 



Qa + aQs 



(8) 



Qi + al —al 
—al Q2 + al 

In the language of physics, a represents the coupling constant of the interaction between 
the networks. Consistently with the rest of the paper, this system will also be referred to 
as the mean-field model of the diagonal interlinking strategy. Regardless of its origin, this 
system exhibits some interesting properties worth discussing. 

Let ^Ni,^Ni-i, •••, ^1 be the set of eigenvectors for the Laplacian of the single network Ai, 
and uni,ujni-i, ■■■,^1 be their relative eigenvalues. Since the perturbation Qb commutes 
with Qa, all the eigenvectors of the interdependent graph are kept unchanged 26|]. All 
the (unperturbed) eigenvalues are degenerate in pairs and, hence, one may define a set of 
eigenvectors based on those of the single networks: 



6 



^2i = 






^21+1 - 


= 




C. 

6._ 



(9) 



The eigenvalues for the total non-interacting system (i.e. a = 0) are the same as for the 
unperturbed system /i2i = fi2i+i = ^i, hence, the ascending sequence of eigenvalues for 
the non-interactive system is u^-^ = 0,0,U]\fj^^i,UNj^^i . . . ,Ui,Ui. When the interaction is 
switched on (i.e. a ^ 0), the even eigenvalues are kept unaltered, while the odd ones increase 
linearly with 2a, 



f^2i = ^ii 

l^2i+i =uJi + 2a. 
For a close to zero, the eigenvector ranking is kept unchanged ^n = <^Ni = 0,fiN-i 



(10) 



C(, fJ'N-2 = ^Ni-1, fJ'N-3 = ^Ni-1 + 2a, . . . , oji, Ui + 2a. However, when a > 



f^iVi -1 



the second 



and third eigenvalues of the interdependent network {^n~i and ^n~2) swap. Therefore, the 
first non-zero eigenvalue increases linearly with 2a up to the value of the isolated networks 
ojn-i-i at which it reaches a plateau. In other words, when a is greater than the threshold 
ai = —3^ the interactive system is capable of synchronizing with the same swiftness as 
the single isolated network. Thus when the system intercommunication channel is quicker 
than the proper time (the inverse of the algebraic connectivity), then the proper time of 
the interactive system equals that of the single network. The critical value aj for the exact 
model corresponds to a critical value of links // to be included to achieve the swiftness of 
the single network: 



// = ajNi 



un~i ■ Ni 



(11) 



If we interpret network robustness as the ability of a system to perform its function upon 
damage or attacks, then it is worth discussing what happens when two networks, Ai and 
A2, originally fully connected by diagonal interlinking B, are subject to some interlink loss. 
Our simple, exact model shows that when these two fully connected networks are subject 
to minor interlink loss, the response of the total interacting system A + aB takes place at 
the same speed as the single component network Ai. In other words, when the operability 



of the control channel via a is mildly reduced, the global system synchronizability does not 
decrease. However if the operability of the connection devices degrades below the critical 
value «/, the synchronization process starts to slow down. From the mean-field approach 
point of view, this means that the system may lose a fraction of interlinks while keeping its 
synchronization time unchanged. 

Following the statistical variant, the parameter a can be regarded as a coupling constant 
or inverse temperature. If one identifies the Fiedler eigenvalue /iAr_i with the internal energy 
of a thermodynamical system, then its first derivative exhibits a jump from zero to a finite 
value. Nevertheless, this derivative does not diverge as expected for a second order transition 



28| 



On the other hand, if one employs the Fiedler eigenvalue as a metric for the synchroniz- 
ability and regards it as a thermodynamical potential such as the free enthalpy, its Legendre 
transform corresponds to the internal energy and exhibits a discontinuity at a = a/. In this 
perspective, one may interpret the observed abrupt change as a first order phase transition. 
Despite this interesting parallel, it is worth noting that the Fiedler eigenvalue and its Legen- 
dre transform are not extensive quantities and, hence, they cannot be properly regarded as 
thermodynamical potentials. However, the behavior of the system closely resembles a phase 
transition. 

To understand the intimate nature of the phase transition, one may inspect the topological 
properties of the eigenvectors. Below the critical value a/, the cut links associated to the 
Fiedler partition lay outside the originally isolated networks (i.e. interlinks are cut), whereas 
just above the critical value, all cut links lay inside the originally isolated networks (i.e. 
intralinks are cut). This means that, below a/, the synchronization is dominated by the 
intralinks in aB^ while beyond «/ the synchronization involves the whole system, A + aB. 

B. General interlinking 

A second important example that may be treated algebraically corresponds to the mean- 
field approximation of the general interlinking strategy. The mean-field approach consists 
of studying a graph with two identical networks interacting via A'^^ weighted connections. 
The interdependence matrix is a matrix with all unitary components: B\2 = J, where J is 
the all ones matrix; the weight of each interlink is a = j^, and 

8 



FIG. 1. Two graphs with 6 black and white nodes respectively and 7 links each are progressively 
interconnected with a) 1 interlink, b) 2 interlinks and c) 3 interlinks. Adding 1 or 2 interlinks 
causes the Fiedler eigenvector to split (depicted by the rectangles) the network into the natural 
partitions Gi and G2 ■ For both cases, the confining links match the added interlinks (dashed 
lines). However, when adding 3 interlinks the Fiedler partition experiences a brusque shift, causing 
the intralinks of the single networks to become the confining links. Thus the added interlinks 
become a part of the Fiedler partitions. 



Q = Qa + aQs 



(12) 



Qi + aNiI -aJ 
-aJ Q2 + aNiI 

As in the previous case, the Qb matrix commutes with Qa and hence a common set of 
eigenvectors can be chosen as in ([9]). The null eigenvalue /iat is always present, while all 
the others experience some increase for a non-trivial a: all eigenvalues /Xj for i smaller than 
A^ — 1, increase for a fixed amount aNi, while /xtv-i increases by twice that quantity. 



N 


= 0, 




N 


_i = 2aiVi, 




i 


= uJi + aNi, ioT i < N - 


-1. 



(13) 



This different rate of growth again implies that there exists a critical value aj beyond 
which the second and third eigenvectors {fJ^N-i and ^n-t) swap. The threshold aj can be 
easily calculated imposing the crossing condition ixn-i = I^n-2'- 



aj 



IvT' 



(14) 



With a = j^, the critical number of links for the general interlinking strategy can be 
also estimated in the mean-field approximation: 



I J = ajNl = ujn-1 ■ Ni 



(15) 



It is worth noting that, the critical number of interhnks corresponding to the mean- field 
theory of the diagonal flTTj) . and general flTSj) interlink strategies, differ simply by a factor of 
2. 



IV. APPROXIMATING hn^i USING PERTURBATION THEORY 

The problem consists in finding the minimum of the associated quadratic form in the 
unitary sphere [x'^x = 1), with the constraint u^x = 0. 



/i = IJ,N-1 



inf 



x'^Qx 



x^O,uTx=0 x'^X 



(16) 



In our case, the matrix Q is the sum of a matrix Qa linking only nodes inside the same 
net, and a "perturbation" aQs that only connects nodes in different networks (Qa + o^Qb)- 
Therefore, we want to find the minimum that satisfies the spectral equations: 



{Qa + qQb - fJ'I)x = 0, 



x'^x = 1, 



u^x = 0. 



(17) 



When the solution is analytical in a, one may express fi and x by Taylor expansion as 



fc=0 



X = y x^^^'a^ 



fe=0 



(18) 
(19) 



Substituting the expansion in the eigenvalue equation f lT6|) gives the hierarchy of equations: 






'X"- 



^J'^o^^'^-^^xW for all k, 
fork> 1, 

for all k. 



(20) 
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A. Explicit approximations up to the second order 



The zero order expansion just provides a simple set of equations: 



a;(o)^(o) ^ i^ (21) 

M^x(o) =0. 

Let {fiN~i)jii , {fiN-i)A2 ^^^ {xn~i)ai 5 {^n~i)a2 denote the smallest non-zero eigenvalue 
and the corresponding eigenvector of Qi, Q2, respectively. Similarly 



(22) 



(x^vJai = l/v^(l,l,...,l,0,0,...,0), 

{xn,)a2 = l/ViV2(0,0,...,0,l,l,...,l). 
will represent the null eigenvectors of network Gi and G2, respectively. When the networks 
are put together, any combination of the former is a null eigenvector. Two special combi- 
nations are worth employing: the trivial solution corresponding to the constant vector: 



1 



Xn 



(1, 



r 




Ni 



i^N^. 



Al 




{XN2) 



A2- 



(23) 



^N" ' ' ' \ N' "''"" V N 
and the other combination orthogonal to the former that represents a useful starting 
point for the perturbation theory: 



X 



(0) 
Af-l 



X 



(0) 



N 



1,-1, 



-1) 




-^(a^TvjAi 



!N2, 
-JT\XN.2)A2- 



(24) 



which satisfies the zero order approximation (!2T]) . 
Fiedler eigenvalue is then null: 



The zero order approximation to the 



/i(o) = 0. 
The first order approximation equations follow from fl2Ul) as: 



(25) 






(^(0))- ^, 



^(1) = (26) 

M^x(^) = 0. 

Taking the projection over x^°^ of the first equation of (126|) . one obtains the first order 
correction /x*^^-* that depends on the zero order eigenvector only: 
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^(i) = (a;(o))^aQ^a;(°) (27) 

A simple case to analyze is that where only one interlink joins Ai with A2: (-Bi2)jj = ^ik^kj', 
in this case {di)kk = ^ik and (^2)^ = 5ji and the perturbation estimate gives: 

/^^^^ = (^(1 + 1) + 1)(^.)' = 1^ > ^^N^l{Q)■ (28) 

where 77 is the single net (A'^i dimensional) unitary vector 77 = l/viVi(l, 1, . . . , 1). When 
k interlinks are included, Qb is just the sum of k contributions of the previous type thus 
^(1) _ 2k_^ That is, the first order correction to the Fiedler eigenvalue increases linearly with 
the number of interlinks. The first order correction to the eigenvector can be evaluated from 
f l26|) as a solution of the linear equation: 

Q^x« = -(agB-/i«)x(°). (29) 

where the operator Qa is invertible out of its kernel {Qav = 0); since [aQs — /i^^^) x^'^^ is 
orthogonal to the kernel, (l29l) is solvable. 

The second order equations follow from (!20|) as 

gAx(2) + aQBX^''^ = /iWx(2) + ;x(i)x(i) + /i(2)xW 
(a;{0))^^(2) + (2;(i))^xW + (a;(2))^xW = (30) 

MX(2) = 

that is, the second order correction is quadratic and equals: 

/x(2) = {x^'^faQn {x^'^) = - {x^'^f Qa {x^'^) < 0. (31) 

As expected /i^^-* is negative, thus improving the estimate of the algebraic connectivity. The 
former perturbation estimates are illustrated in Fig. [2] together with numerical simulations. 
Perturbation theory may also be applied to any initial eigenvector of the unperturbed 
networks. Different perturbations aB will have different effects on the quadratic form of 
( IT6II associated with all initial eigenvectors. Therefore, it may happen that the perturbed 
value of /i obtained starting from x^^^ is smaller than the quadratic form associated with the 
xn^i (the unperturbed eigenvector in (|9])) or some other educated guess. This is precisely 
the origin of the phase transition. 
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RR p.theory n +^i 
BAp.theory n'^'V 
WS p.theory fi'^'+^ 
LAp.theory n V 

- RR simulations 

- BA simulations 

- WS simulations 
LA simulations 
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^ 




interlinks interlinl<s 

(a) Diagonal strategy (b) General strategy 

FIG. 2. Simulated (solid lines) and theorized (dashed lines) algebraic connectivity fiN-i of four 
graph models with N = 1000 nodes, as interlinks are added between single networks following 
two strategies: diagonal interlinks (left image) and general interlinks (right image). Perturbation 
theory best approaches fJ-N-i for the parabolic region of the diagonal interlinks strategy, which 
saturates after adding "^^' links, as we detailed in section illli 

The estimates resulting form the second order perturbation theory are compared in Fig. [2] 
with the results of numerical calculations. As can be seen, for both the diagonal and the 
general strategies the agreement is good up to the phase transition where the starting point 
of the perturbation theory should be changed. 

B. Perturbative approximations and upper bounds 

Since we are dealing with a constraint optimization problem, finding a minimum of a 
positive form, any test vector v provides an upper bound for the actual minimum value: 



/i = fiN-i < r 



V^ V 



(32) 



The perturbation theory provides natural candidates as test vectors. The zero order 
solution provides the simplest inequality: 



fJ'N^iiQ) < a 



^3,(0))Tg^(0) 



an 



(1) 



(x(o))^x(o) 
The first order approximation provides a better (i.e. lower) upper bound: 



(33) 
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that is: 



/^^-i(Q) < 1 + aHxWY • (^^^ 



(1) 



which for small enough a is always lower than afi 



V. SIMULATIONS 

Previous sections provided basic means to understand the dependence of the algebraic 
connectivity on the topology of the interdependence links. In this section we will intro- 
duce model networks to test the predictability and the limits of the mean-field and the 
perturbation approximations. 

A. Interdependent networks model 

Our interdependent network model consists of two main components: a network model 
for the single networks, and the rules by which the two networks are linked. In other words, 
to model two interdependent networks one needs to select two model networks and one 
interlinking strategy. 

In the numerical simulations discussed here, we considered four different graph models 
for our coupled networks. These models exhibit a wide variety of topological features and 
represent the four different building blocks: 



Random Regular (RR): random configuration model introduced by Bollobas [29|. 
All nodes are initially assigned a fixed degree di = k,i & N . The k degree stubs are 
then randomly interconnected while avoiding self-loops and multiple links. 

Barabasi-Albert (BA): growth model proposed by Barabasi et. al. [l| whereby new 
nodes are attached to m already existing nodes in a preferential attachment fashion. 
For large enough values of A^, this method ensures the emergence of power-law behavior 
observed in many real-world networks. 
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• Watts-Strogatz (WS): randomized circular lattice proposed by Watts et. al. |2| 
where all nodes start with a fixed degree k and are connected to their | immediate 
neighbors. In a second stage, all existing links are rewired with a small probability 
p, which produces graphs with low average hopcount yet high clustering coefficient, 
which mimics the small-world property found in real-world networks. 

• Lattice (LA): a deterministic three-dimensional grid which loops around its bound- 
aries (i.e. a geometrical torus). 

The input parameters for each model are set such that all graphs have the same number of 
nodes and links. In addition, all simulated graphs consist of a single connected component, 
i.e. random graphs containing more than one connected component were discarded. 

We define two strategies to generate the interdependency matrix 5, which we analytically 
solved in section lllli 

• diagonal interlinking strategy: links are randomly added to the diagonal elements of 
-B, thus linking single network's analogous nodes. 

• general interlinking strategy: random links are added to B without restrictions, gen- 
erating a random interconnection pattern. 

In the next sections, we explore the effects of the two interlinking strategies on the Fiedler 
partition. However, some results can be extended to more complicated situations, including 
different synthetic networks or more complex linking strategies. 

B. Partition quality metrics 

Let us introduce some preliminary definitions, required for understanding of our numerical 
results. We define a graph bipartition of G as the two disjoint sets of nodes {7^,5}, where 
TZVJ S = M . We define the natural partition of G as the partition with the two original 
node sets: TZ = A/i, and S = A/2. The number of nodes in 7?. and S is counted by their 
cardinality \TZ\ and |iS|, respectively. In addition, we express the number of links with one 
end node in 7?. an another end node in 5 as / {TZ, S) = I {S, TZ). Fiedler partitioning bisects 
the nodes in M into two clusters, such that two nodes i and j belong to the same cluster 
if XiXj > 0, i.e. the corresponding components of the Fiedler eigenvector x have the same 
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FIG. 3. The four main partition sets are displayed: A/i (set of black nodes), ^2 (set of white nodes), 
TZ (set of nodes within the gray rectangle), and S (set of nodes within the white rectangle). Both, 
the partition sets and the interlinks (dashed lines) were arbitrary chosen for illustration purposes 
and do not represent the corresponding Fiedler partition. 

sign. For example, if the coupling strength a in ([3]) is zero, the bipartition resulting from 
Fiedler partitioning is equivalent to the two natural clusters, i.e. TZ = Gi and S = G2- 

The intersection between the two single graphs (A/i, A/2) and the Fiedler partitions {IZ, S) 
of the interdependent network yields four node subsets, defined as follows and illustrated in 
Fig. [HI (a) TZi as the intersection between the positive Fielder partition and A/i, (b) Si as 
the intersection between the negative Fielder partition and A/i, (c) 7^2 as the intersection 
between the positive Fielder partition and A/2, (d) S2 as the intersection between the positive 
Fielder partition and N'2- By construction the four defined groups are disjoint, and the union 
of all groups equals the full set of nodes. 

In order to quantify the properties of the Fiedler partition, we study the following set of 
metrics: 

• Fiedler cut-size = ^ \ . It represents the fraction of links with one end in 7?. and 
another end in S (irrespective of the directionality of the link) over the starting number 
of links. 



interdependence angle, defined as the angle between the normalized Fielder vector x 
and the versor x^°\ introduced in (!23l) . The interdependence angle is minimized when 
the Fiedler vector is parallel to the natural partition, i.e. x^^\ 
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entropy of the squared Fiedler vector components = — ^j=i Xj^ logXj^. Based on 
Shannon's information theory metric, the entropy indicates how homogeneous the 
values in x are, similarly to the participation ratio or vector localization. The higher 
the entropy, the lower is the dispersion among the values in x. 

Some partition quality metrics may be undefined if the Laplacian matrix Q is defective 



30[. In particular, if the second and third largest eigenvalues of Qa + aQs are equal 
^N-i = f^N~2 then any linear combination x' = axjy^i + &X7V-2 is also an eigenvector of 
Q with eigenvalue fiN-i, thus the Fiedler vector is not uniquely defined. However we will 
ignore these cases, which tend to occur only in graphs with deterministic structures (e.g. 



the cycle graph |26|). 



C. Diagonal Interlinking Strategy 

1. Strategy Description 

The diagonal interlinking strategy consists of adding links between the respective compo- 
nents of two identical networks. We can add as little as 1 link and as many as A^ links. This 
strategy was chosen to achieve the maximum effect by meticulously adding a small number 
of interlinks. A simplified physical example would be that of two flat metal plates: a hot 
one, and a cold one. If the objective is to equalize their temperature as fast as possible, we 
should adjust the plates side by side so as to maximize the heat transfer, which is equivalent 
to the diagonal interlinks strategy. 

2. Initial and final states 

We will refer to the natural, initial or unperturbed state as the scenario where there exist 
no interlink connecting the two networks Gi and G2- The left hand side of Fig. [1] and Fig. H] 
illustrate the network configuration and partition quality metrics, respectively. In this case 
the algebraic connectivity dips to a null value, as no communication is possible; the Fiedler 
partition then becomes undetermined. However, the sign of x*-*^^ in (l29l) allows splitting the 
network into two clusters V = Gi and Q = G2, corresponding to the isolated component 
networks. 
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FIG. 4. (Color online) Four metrics' averages are displayed to evaluate the effect of adding interlinks 
following the diagonal strategy: algebraic connectivity (hn-i), Fiedler cut {l(TZ,S)/Li + L2), 
interdependence angle (acos (x x^'^Y ||x|| ||x'''^||)), and entropy (— X]j=i Xj^ logXj^). All metrics 
experience a transition that sharpens for increasing N. BA and RR graphs transition around 80% 
added interlinks, whereas WS and LA graphs transition around 20%. The size of the network 
A^i has a relatively little impact on BA and RR curves, which suggests that the transition is 
independent of the network size A'^i. The flat lines signaled with arrows in the top left plot 
benchmark the average algebraic connectivity of the A'^i = 1,000 respective single networks. 



The final state of the diagonal interlink strategy corresponds to 100% or A^ added in- 
terlinks, thus B = I and the Fiedler vector becomes the vector {x]\f^i{Qi),X]\f^i{Qi)), as 
demonstrated in section UTTI The final partition depends exclusively on Gi and G2, indepen- 
dently of B. Since we assume Gi = G2, the final cut consists purely of a subset of intralinks 
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FIG. 5. (Color online) The fluctuation a"^ of the Fiedler cut and the interdependence angle are 
displayed to evaluate the effect of adding interlinks following the diagonal strategy. The narrowing 
peaks illustrate the sharpening of the transition observed in Fig. HI 

of Ai, as illustrated in the right hand side of Fig. [1] By adding (or removing) interlinks 
between the two independent networks we switched from a purely interlink cut to a purely 
intralink cut, which is the essence of the phase transition. 

The sudden partition quality metric shifts reflect the change observed in the Fiedler 
partition, while providing evidence for the existence of the transition. This is testified by 
the narrowing of the region where the shifts occurs while increasing (or decreasing) A^ in 
Fig. m Similarly, the fluctuations of the same quantities exhibit shrinking peaks as can be 
seen in Fig. [51 Upon increasing the size A^ of the system, the transition point seems to 
approach an asymptotic value. As discussed in section llllt the mean-field theory predicts 
this critical value to be // = ^^-iWiJ- i _ gjnce jJiN-i{Qi) has a non-trivial lower bound 



for increasing A^ 31[, and its algebraic connectivity is kept unchanged by the perturbation, 
there will be a critical number of links beyond which fiN-iiQ) does not change. This critical 
point corresponds to the transition from an interlink cut to an intralink cut. 

The precise location of the jump in the simulated experiment, i.e. the critical value of 
interlinks per node, depends on the graph model. However, the phase transition is a general 
phenomenon, which occurrence only depends on the fact that there exists a Fiedler cut for 
the single networks. 
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3. Effect on partition quality metrics 

We further investigate the properties of the phase transition by looking at how partition 
metrics in Fig. H] evolve as interlinks are added to the B matrix. 

The algebraic connectivity starts at its minimum value ~ ^^ as predicted by (!28|l . 
which grows until it reaches its maximum value fiN-i{Qi) when sufficient interlinks are 
added. This means that a network with 100% diagonal interlinks and the same network 
with 90% interlinks synchronize virtually at the same speed. Comparing the final values of 
the algebraic connectivity, it is remarkable that random networks synchronize faster than 
lattice networks. This is reasonably due to the longer average distance in the latter. 

The Fiedler cut starts at ^ for a single added interlink. Notice that it increases linearly 
with the percentage of interlinks, because all added interlinks directly become part of the 
Fiedler cut. For all networks, we observe a tipping point (which depends on the network 
type) upon which adding a single link abruptly readjusts the partition: the Fiedler cut 
switches from pure interlink cutting to a cutting of an invariable set of intralinks. This 
abrupt change breaks the linearity. 

The interdependence angle metric tells us that the Fiedler vector starts being parallel 
to the first order approximation x^'^^ for 1 added interlink. Progressively, the Fiedler vector 
crawls the A^-dimensional space up to the transition point, where it abruptly jumps to 
the final (orthogonal) state {xiy^i{Qi) , xn-i{Qi)) ■ Similarly to the interdependence angle, 
the high values of entropy reflect the flatness of x^^\ where all components have (almost) 
the same absolute value. At this initial point, entropy is maximum and almost equal to 
log{2Ni), which tells us that the initial partition consists purely of interlinks. When the 
partition turns to the final state, the entropy is instantly shaped by the network topologies 
of Ai thus dropping to relatively much lower values. Notice that, for all values of A^, the 
highest final entropy is attained by the lattice graph due to its regular structure, as seen in 
Fig.H 

4- Network Model Differences 

RR and BA synchronize relatively faster than deterministic networks because random 
interconnections shorten the average hopcount, thus bringing all elements of the network 
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closer (creating the small-world effect |2]). For the particular case of BA, we observe the 
emergence of a dominant partition which contains approximately 90% of the total number 
of nodes. 

There exists a significant difference between the 1,000 node lattice and the 10,000 node 
lattice, which is expected due to the variable size response of network models. We conjecture 
that this difference is caused by the average geodesic distance: the average node distance 
for a three dimensional lattice lattice graphs grows with v^, as opposed to random models, 
which usually display logarithmic increases. Interestingly, Fig. |4] shows that small lattices 
synchronize faster than WS, but the situation is soon reversed for higher A^. 

To test whether the phase transition is merely an artifact of our synthetic models, addi- 
tional simulations were carried out using real topologies from the KONECT dataset. Simu- 
lations verify that the transition from the natural partition to the final orthogonal partition 
also occurs in real networks. However, the transition takes place very early in the link 
addition process, due to the poor synchronization capabilities of networks not designed for 
such purpose. The interpretation of such result is that, to provide that real network with 
a complete backup mirror without synchronization delays, a small number of interlinks are 
required. 

D. General Interlinks Strategy 

1. Strategy Description 

As a variation of the localized diagonal interlinking strategy, our second strategy ran- 
domly draws interlinks among any pair of nodes belonging to different networks. Mean-field 
approximation provides us with exact results, however perturbation analysis loses its power 
when too many links are added, i.e. the perturbation can no longer be regarded as small. 
For this reason, we cannot predict an exact asymptotic state as for the diagonal strategy. 
We have limited our simulations to the inclusion of up to 4 interlinks per node. 

2. Effect on partition quality metrics 

We can observe that the algebraic connectivity of all models experiences two regimes, 
upon the progressive addition of interlinks as illustrated in Fig. |6l Initially, for a small 
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FIG. 6. (Color online) Four metrics' averages are displayed to qualitatively evaluate the 
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transitions are not as sharp as in the diagonal strategy scenario. The flat lines signaled with arrows 
in the top left plot represent the average algebraic connectivity of the A'^i = 1, 000 respective single 
networks. 



number of added links, the initial state dips to a minimum as is the case for the diagonal 
strategy and represents a good starting point for the perturbation theory. As we increase the 
number of interlinks, the average algebraic connectivity and Fiedler cut curves show a linear 
increase. At the critical number of links Ij = fiN-i ■ N, the average slope switches regime 
by damping to half its value, as seen in Figj6a] and Fig. l6b| which is in perfect agreement 
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FIG. 7. (Color online) The fluctuation a^ of the Fiedler cut and the interdependence angle are 
displayed to evaluate the effect of adding interlinks following the general strategy. The increasing 
Fiedler cut's fluctuations do not hint the existence of a transition. However the peaking fluctuation 
displayed by the interdependence angle suggests the presence of a narrowing transition. 



with our theoretical prediction ( TTSl) . However not only the average, but also the fluctuations 
steadily increase, as illustrated by Fig. [Tal High fluctuations are expected, due to the large 
set of available graph configurations. 

As we can see from the interdependence angle in Fig. [6l in the first regime the natural 
partition is partially preserved up to Ij. The interdependence angle experiences a sharp 
increase at the turning point, which further narrows as A^ increases as seen in Fig. [70 This 
is due to the fact that the Fiedler cut in all our isolated model networks scales less than 
linearly with the network size, which is consistent with the picture of a phase transition 
between a Fiedler cut dominated by interlinks and an other dominated by intralinks. As 
opposed to the diagonal strategy, the final eigenvector is not strictly identical to the Fiedler 
eigenvector of the isolated networks xjv_i, but it also involves interlink cuts. This is due to 
the fact that in the general case x^-\ does not belong to the kernel of Qb as opposed to the 
diagonal case. 

The exact location of the phase transition can also be predicted employing perturbation 
theory, by imposing the perturbed value Xfq-\{ff) of the configuration to be equal to that 
achieved starting from the Xfq^\[Q\) initial state. However, the resulting formulas are not 
particularly simple and their numerical calculation requires a time comparable with the 
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Fiedler eigenvalue evaluation of the sparse metrics. For this reasons such estimates are not 
reported here. 

3. Network Model Differences 

Let us focus on the case of adding a small number of interlinks in the range [1, iV]. The 
diagonal strategy will synchronize faster than the general strategy in the case of RR and 
BA, as illustrated in Fig. [6al On the other hand, the general strategy synchronizes faster 
in WS and LA models. Thus if we were to add precisely A'^i general interlinks between two 
identical networks, regular structures would (relatively) benefit the most. 

For BA, the fraction of intralinks belonging to the Fiedler partition decreases with in- 
creasing number of interlinks, whereas the ^ ratio increases. This hints that nodes group 
into high degree clusters (with a high link/node ratio) and a low degree clusters (with a low 
link/node ratio). In addition, BA's entropy experiences the highest drop, which indicates 
that the Fiedler vector is highly localized around a small set of nodes. 

The difference between random and grid networks still exists for the general strategy, 
but it is not as predominant as in the diagonal case. This effect is expected due to the 
randomization resulting from the random addition of links to regular structures, which is 
the conceptual basis of the WS model. In general, we observe that the optimal link addition 
strategy depends on the network topology. 

VI. CONCLUSIONS 

This paper aims to provide general results concerning the synchronization of interdepen- 
dent identical networks. We provided evidence that upon increasing the number of interlinks 
between two originally isolated networks, their synchronizability experiences a phase transi- 
tion. That is, there exists a critical number of diagonal interlinks beyond which any further 
inclusion does not enhance synchronization capabilities at all. Similarly, there exists a criti- 
cal number of general interlinks beyond which algebraic connectivity increments at half the 
original rate. 

The exact location of the transition depends exclusively on the algebraic connectivity of 
the graph models, and it is always observed regardless of the interconnected graphs. For the 
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two proposed interconnection strategies, the critical number of interlinks that triggers the 
transitions can be predicted correctly by mean-field approximations : ^^-i'-^^-'' ^ links for 
the diagonal interlinks strategy, and fiN-iiQi) ■ Ni links for the general interlinks strategy. 
By resorting to perturbation theory we have provided upper bounds for the total algebraic 
connectivity of the interdependent system and means to estimate it. 

This paper beacons a significant starting point to the understanding of the mutual net- 
works synchronization phenomena, as we have just started studying this extremely interest- 
ing field. Nonetheless, different linking strategies should be researched and general theory 
developed. Regarding the mutual synchronization of heterogeneous networks (i.e. Ai ^ A2), 
preliminary results confirm the existence of phase transitions with similar features to the 
general random linkage of identical networks. However, we could not observe any dominant 
strategy as in the case with the diagonal interlinking. 
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